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Abstract. A new fluid-dynamic model is developed to nu¬ 
merically simulate the non-equilibrium dynamics of polydis- 
perse gas-particle mixtures forming volcanic plumes. Start¬ 
ing from the three-dimensional N-phase Eulerian transport 
equations ( |Neri et ah 2003| l for a mixture of gases and 
solid dispersed particles, we adopt an asymptotic expansion 
strategy to derive a compressible version of the first-order 
non-equilibrium model (Ferry and Balachandar 200 l| l, valid 
for low concentration regimes (particle volume fraction less 
than 10“^) and particles Stokes number (St, i.e., the ratio 
between their relaxation time and flow characteristic time) 
not exceeding about 0.2. The new model, which is called 
ASHEE (ASH Equilibrium Eulerian), is significantly faster 
than the N-phase Eulerian model while retaining the capa¬ 
bility to describe gas-particle non-equilibrium effects. Direct 
numerical simulation accurately reproduce the dynamics of 
isotropic, compressible turbulence in subsonic regime. For 
gas-particle mixtures, it describes the main features of den¬ 
sity fluctuations and the preferential concentration and clus¬ 
tering of particles by turbulence, thus verifying the model re¬ 
liability and suitability for the numerical simulation of high- 
Reynolds number and high-temperature regimes in presence 
of a dispersed phase. On the other hand, Large-Eddy Nu¬ 
merical Simulations of forced plumes are able to reproduce 
their observed averaged and instantaneous flow properties. In 
particular, the self-similar Gaussian radial profile and the de¬ 
velopment of large-scale coherent structures are reproduced, 
including the rate of turbulent mixing and entrainment of 
atmospheric air. Application to the Large-Eddy Simulation 
of the injection of the eruptive mixture in a stratified atmo¬ 
sphere describes some of important features of turbulent vol¬ 
canic plumes, including air entrainment, buoyancy reversal, 
and maximum plume height. For very fine particles (St —0, 
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when non-equilibrium effects are negligible) the model re¬ 
duces to the so-called dusty-gas model. However, coarse 
particles partially decouple from the gas phase within ed¬ 
dies (thus modifying the turbulent structure) and preferen¬ 
tially concentrate at the eddy periphery, eventually being lost 
from the plume margins due to the concurrent effect of grav¬ 
ity. By these mechanisms, gas-particle non-equilibrium pro¬ 
cesses are able to influence the large-scale behavior of vol¬ 
canic plumes. 


1 Introduction 


Explosive volcanic eruptions are characterized by the injec¬ 
tion from a vent into the atmosphere of a mixture of gases, 
liquid droplets and solid particles, at high velocity and tem¬ 
perature. In typical magmatic eruptions, solid particles con¬ 
stitute more than 95% of the erupted mass and are mostly 
produced by the brittle fragmentation of a highly viscous 
magma during its rapid ascent in a narrow conduit ( |Wil-| 
1976t Sparks[ 19781, with particle sizes and densities 


son 


spanning over a wide range, depending on the overall char¬ 
acter and intensity of the eruption ([Kaminski and Jaupart 


1998| Kueppers et al. 2006)1. The order of magnitude of the 


plume mixture volumetric concentration very rarely exceed 
Cs ~ 3* 10“^, because the order of magnitude of the ejected 
fragments density is Ps ^ 10^ kg/m^. Thus, the plume mix¬ 
ture con be considered mainly as a dilute suspension in the 
sense of Elghobashi (1991 |1994| l. This threshold for Cg is 
overcome in the dense layer forming in presence of pyroclas¬ 
tic density currents (see e.g. |Orsucm]|2014| l. In the literature, 
collisions between ash particles are usually disregarded when 
looking at the dynamics of volcanic ash plume, because this 
dilute character of the plume mixture (cf.jMorton et al.||1956 
[Woods) |20T0 | i. 

After injection in the atmosphere, this multiphase eruptive 
mixture can rise convectively in the atmosphere, forming ei- 
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ther a buoyant volcanic plume or collapse catastrophically 
forming pyroclastic density currents. Since these two end- 
members have different spatial and temporal scales and dif¬ 
ferent impacts on the surrounding of a volcano, understand¬ 
ing the dynamics of volcanic columns and the mechanism of 
this bifurcation is one of the topical aims of volcanology and 
one of the main motivations for this work. 

The term volcanic column will be adopted in this pa¬ 
per to genetically indicate the eruptive character (e.g. con¬ 
vective/collapsing column). Following the fluid-dynamic 
nomenclature, we will term jet the inertial regime of the 
volcanic column and plume the buoyancy-driven regime. A 
forced plume is characterized by an initial momentum-driven 
jet stage, transitioning into a plume. 

In this work, we present a new computational fluid- 
dynamic model to simulate turbulent gas-particle forced 
plumes in the atmosphere. Although the focus of the paper is 
on multiphase turbulence in subsonic regimes, the model is 
also suited to be applied to transonic and supersonic flows. 
In many cases, indeed, the eruptive mixture is injected in 
the atmosphere at pressure higher than atmospheric, so that 
the flow is initially driven by a rapid, transonic decompres¬ 
sion stage. This is suggested by numerical models predicting 
choked flow conditions at the volcanic vent ( |Wilson| |1980t 
Wilson et al. 1980|l, implying a supersonic transition above 


the vent or in the crater ( |Kieffer| |1984 [ [Woods and Bower| 
1995 1 Koyaguchi et al. 2010|l and it is supported by held 


evidences of the emission of shock waves during the initial 
stages of an eruptions ( |Morrisse3^ 19971. Despite the impor¬ 
tance of the decompression stage for the subsequent devel¬ 


opment of the volcanic plume (Pelanti and LeVeque 2006 


et al. 


Ogden et al. 2008b[ Orescanin et al. 2010 Carcano et al. 


2013 I and for the stability of the eruptive column (Ogden 
2008a|l, our analysis is limited to the plume region 


where flow pressure is equilibrated to the atmospheric pres¬ 
sure. From laboratory experiments, this is expected to occur 
within less than 20 inlet diameters above the ground ( jYuceilj 
and Otiigen 20021. 


1.1 Dusty gas modeling of volcanic plumes 

Starting from the assumption that the large-scale behavior 
of volcanic columns is controlled by the bulk properties of 
the eruptive mixture, most of the previous models of vol¬ 
canic plumes have considered the eruptive mixture as homo¬ 
geneous (i.e., they assume that particles are perfectly coupled 
to the gas phase). Under such hypothesis, the multiphase 
transport equations can be largely simplified and reduce to a 
set of mass, momentum and energy balance equations for a 
single fluid (named dusty-gas or pseudo-gas) having average 
thermo-fluid dynamic properties (mixture density, velocity 
and temperature) and equation of states accounting for the 
incompressibility of the particulate phase and gas covolume 
( |Marble|[T970l ). 


By adopting the dusty gas approximation, volcanic plumes 

|T96^ 


have been studied in the framework of jet (Prandtl 


and plume theory ([Morton et al.j 1956 Morton 19591. One 


dimensional, steady-state pseudo-gas models of volcanic 
plumes have thus had a formidable role in volcanology to 
identify the main processes controlling their dynamics and 
scaling properties (|Wilson[[1976[[Woods[[1988[[Sparks et al. 
[T997l l. 

Accordingly, volcanic plume dynamics is schematically 
subdivided into two main stages. The lower, jet phase is 
driven by the initial flow momentum. Mixture buoyancy is 
initially negative (the bulk density is larger than atmospheric) 
but the mixture progressively expands adiabatically thanks to 
atmospheric air entrainment and heating, eventually under¬ 
going a buoyancy reversal. When buoyancy reversal does not 
occur, partial or total collapse of the jet from its maximum 
thrust height (where the jet has lost all its initial momentum) 
and generation of pyroclastic density currents are expected. 

Above the jet thrust region, the rise of volcanic plumes is 
driven by buoyancy and it is controlled by turbulent mixing 
until, in the stratified atmosphere, a level of neutral buoyancy 
is reached. Above that height, the plume starts to spread out 
achieving its maximum height and forming an umbrella ash 
cloud, dispersing in the atmosphere and slowly falling-out. 

In one-dimensional, time-averaged models, entrainment of 
atmospheric air is described by one empirical coefficient (the 
entrainment coefficient) relating the influx of atmospheric air 
to the local, vertical plume velocity. The entrainment coeffi¬ 


cient also determines the plume shape (Ishimine 20061 and 


can be empirically determined by means of direct field ob¬ 
servations or ad-hoc laboratory measurements. 

Further development of volcanic plume models have in¬ 
cluded the influence of atmospheric stratification and humid- 


ity (Woods |1993 

Glaze and Baloga| 1996|l, the effect of 

cross wind ([Bursi 

c 200 Ijl, loss and reentrainment of solid 

particles from plume margins (IWoods and Bursik 

1991 

Veitch and Woods 

2002|), and transient effects ( Scase[ 

2009 


Woodhouse et al.[ 2015|l. However, one-dimensional models 


strongly rely on the self-similarity hypothesis, whose validity 
cannot be experimentally ascertained for volcanic eruptions. 

To overcome the limitations of one-dimensional models, 
three-dimensional dusty-gas models have been developed to 


simulate volcanic plumes. Suzuki (2005i have developed 
a three-dimensional dusty gas model (SK-3D) able to ac¬ 
curately resolve the relevant turbulent scales of a volcanic 
plume, allowing a first, theoretical determination of the en¬ 
trainment coefficient ( [Suzuki and Koyaguchi[|20i0) l, without 
the need of an empirical calibration. 

To simulate the three-dimensional large-scale dynamics of 
volcanic plumes including particle settling and the complex 
microphysics of water in volcanic plumes, the ATHAM (Ac¬ 
tive Tracer High Resolution Atmospheric Model) code has 
been designed (Oberhuber et al.[[l998[[Graf et al.[[T999[[Van| 


Eaton et al. 20151. ATHAM describes the dynamics of gas- 
particle mixtures by assuming that particles are in kinetic 
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equilibrium with the gas phase only in the horizontal compo¬ 
nent, whereas along the vertical direction they are allowed to 
have a differential velocity. Thermal equilibrium is assumed. 
In this sense, ATHAM relaxes the dusty-gas approximation 
(while maintaining its fundamental structure and the same 
momentum transport equations) by describing the settling of 
particles with respect to the gas. 

1.2 Multiphase flow models of volcanic plumes 


metical code development are reported in Section]^ Section 
|4] focuses on verification and validation issues in the context 
of applications to turbulent volcanic plumes. In particular, 
here we discuss: three-dimensional numerical simulations of 
compressible, isotropic turbulence (with and without parti¬ 
cles); experimental scale forced plumes; Sod’s shock tube 
problem. Finally, Section [^presents numerical simulations 
of volcanic plumes and discusses some aspects related to nu¬ 
merical grid resolution in practical cases. 


Notwithstanding all the above advantages, dusty-gas models 
are still limited by the equilibrium assumption, which can be 
questionable at least for the coarsest part of the granulometric 
spectrum in a plume. Turbulence is indeed a non-linear, mul¬ 
tiscale process and the time and space scales of gas-particle 
interaction may be comparable with some relevant turbulent 
scales, thus influencing the large-scale behavior of volcanic 
plumes. 

To model non-equilibrium processes, Eulerian multiphase 
flow models have been developed, which solve the full set 
of mass, momentum, and energy transport equations for a 
mixture of gas and dispersed particles, treated as interpen¬ 
etrating fluids. Valentine and Wohletz ( 1989|l and Dobran 


et al. ( 1993|l; Neri and Dobran ( 1994|l have first analyzed the 


influence of erupting parameters on the column behavior to 
identify: By means of two-dimensional numerical simula¬ 
tions, they individuate a threshold from collapsing to con¬ 


vective columns. Lately, two-dimensional (Di Muro et al. 


|2004| IDartevelle et al. 2004| l and three-dimensional numeri¬ 
cal simulations (Esposti Ongaro et al. 2008) 1 has contributed 
to modify the view of a sharp transition between convect- 
ing and collapsing columns in favor of that of a transitional 
regime, characterized by a progressively increasing fraction 
of mass collapsing. However, previous works could not in¬ 
vestigate in detail the non-equilibrium effects in volcanic 
plumes, mainly because of their averaged description of tur¬ 
bulence: a detailed resolution of the relevant turbulent scales 
in three dimensions would indeed be computationally pro¬ 
hibitive for N-phase systems. 

The main objective of the present work is therefore to de¬ 
velop a new physical model and a fast three-dimensional nu¬ 
merical code able to resolve the spatial and temporal scales of 
the interaction between gas and particles in turbulent regime 
and to describe the kinetic non-equilibrium dynamics and 
their influence on the observable features of volcanic plumes. 
To this aim, a development of the so-called equilibrium- 


Eulerian approach (Ferry and Balachandar 2001) Balachan- 
dar and Eaton 2010|l has been adopted. It is a general¬ 


ization of the dusty-gas model keeping the kinematic non¬ 
equilibrium as a first order correction of the |Marble| ( |1970| l 
model with respect to the Stokes number of the solid parti¬ 
cles/bubbles in the mixture. 

The derivation of the fluid dynamic model describing the 
non-equilibrium gas-particle mixture is described in detail in 
Section]^ The computational solution procedure and the nu- 


2 The multiphase flow model 


To derive an appropriate multiphase flow model to describe 
gas-particle volcanic plumes, we here introduce the non- 
dimensional scaling parameters characterizing gas particle 
and particle particle interactions. 

The drag force between gas and particles introduces in the 
system a time scale Tj, the particle relaxation time, which 
is the time needed to a particle to equilibrate to a change of 
gas velocity. Gas-particle drag is a non-linear function of the 
local flow variables and, in particular, it depends strongly on 
the relative Reynolds number, defined as: 


_ Pgl'^s Ug\ds 


( 1 ) 


here ds is the particle diameter, /5g is the gas density, p, is 
the gas dynamic viscosity coefficient and Mg(s) is the gas 
(solid) phase velocity held. Being pg(s) the gaseous (solid) 
phase density and Cs = V^jV the volumetric concentration 
of the solid phase, it is useful to define the gas bulk den¬ 
sity pg = (1 — es)/5g ~ pg and the solid bulk density ps = CsPs 
(even though in our applications Cs is order 10“^, ps is non- 
negligible since Ps/Pg is of order 10^). 

For an individual point-like particle (i.e., having diameter 
ds much smaller than the scale of the problem under analy¬ 
sis), at Res < 1000, the drag force per volume unity can be 
given by the Stokes’ law: 


/d= —(Ug-Ms), 
'T's, 


( 2 ) 


where 


Pg 18 z/(/)c(Res) 


(3) 


is the characteristic time of particle velocity relaxation with 
respect the gas, ps is the particle density, i/ is the gas kine¬ 
matic viscosity and = l-l-0.15Re° ®®^ is a correction fac¬ 
tor (obtained from the Schiller-Naumann correlation) for fl- 


nite particle Reynolds number (cf IClift et al. 

11978) Bal- 

lachandarl |2009t IBalachandar and Eaton 

2010 

ICerminara 


pressure gradient, the added mass, the Basset history and the 
Saffman terms, because we are considering heavy particles: 


Ps/Pg^ 1 (cf. Ferry and Balachandar 2001 Bagheri et al. 
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|2013| l. Equation Q has a linear dependence on the fluid- 
particle relative velocity only when Res ^ 1, so that — 1 
and the classic Stokes drag expression is recovered. On the 
other hand, if the relative Reynolds number Res grows, non¬ 
linear effects become much more important in Eq. Q. The 
Clift et al. ( 1978|l empirical relationship used in this work has 


been used and tested in a number of papers (e.g., |Balachan- 
dar and Eaton 2010[ Wang and Maxey| 1993[ Bonadonna 


et al.| 20021, and it is equivalent to assuming the following 


gas-particle drag coefficient; 

24 


CD(Res) = ^(l + 0.15Re°®®^). 
Res 


(4) 


[Wang and Maxey (19931 discussed nonlinear effects due to 
this correction on the dynamics of point-like particles falling 
under gravity in an homogeneous and isotropic turbulent sur¬ 
rounding. We recall here the terminal velocity that can be 
found by setting ttg = 0 in Eq. (|^ is; 


/ 4dsPs 


(5) 


As previously pointed out, correction used in Eq. Q is valid 
if Res < 10^, the regime addressed in this work for ash parti¬ 
cles much denser then the surrounding fluid and smaller than 
1mm. As shown by [Balachandar ( |2009[ ), maximum values 
of Res are associated to particle gravitational settling (not to 
turbulence). Using formula and (|^, it is thus possible to 
estimate Res of a falling particle with diameter ds- We obtain 
that Res is always smaller than 10^ for ash particles finer than 
1 mm in air. If regimes with a bigger decoupling needs to be 
explored, more complex empirical correction has to be used 


for (fic ( Neri et ak] 2003t Burger and Wendland 20011. 

The same reasoning can be applied to estimate the ther¬ 
mal relaxation time between gas and particles. In terms of 
the solid phase specific heat capacity and its thermal con¬ 
ductivity fcg, we have; 


tt = 


2 

Nus k„ 12’ 


(6) 


where Nus = Nus(Res,Pr) is the Nusselt number, usually 
function of the relative Reynolds number and of the Prandtl 
number of the carrier fluid ( |Neri et ari|2003| l. In terms of tt, 
the heat exchange between a particle at temperature and 
the surrounding gas at temperature Tg per unit volume is; 


Qr = ^(7;-Tg). 


Tt 


(7) 


Comparing the kinetic and thermal relaxation times we 
get; 

Tt 3 20c C's/r 

Ts 2 Nus kg 

In order to estimate this number, firstly we notice that factor 
20c/Nus tends to 1 if Res —0, and it remains smaller than 


( 8 ) 


~ 2 if Res < 10^ ( Neri et al. 2003 Cerminar^ 2015b I. Then, 
in the case of ash particles in air, we have (in SI units) /i ~ 
10“®, Cs ~ 10^, fcg ~ 10“^. Thus we have that tt/t^ — 1, 
meaning that the thermal equilibrium time is typically of the 
same order of magnitude of the kinematic one. This bound is 
very useful when we write the equilibrium-Eulerian and the 
dusty gas models, because it tells us that the thermal Stokes 
number is of the same order of the kinematic one, at least for 
volcanic ash finer than 1 mm. 

The non-dimensional Stokes number (St) is defined as the 
ratio between the kinetic (or thermal) relaxation time and a 
characteristic time of the flow under investigation Tp, namely 
Sts = TsjTp. The definition of the flow time-scale can be 
problematic for high-Reynolds number flows (typical of vol¬ 
canic plumes), which are characterized by a wide range of 
interacting length- and time-scales, a distinctive feature of 
the turbulent regime. Eor volcanic plumes, the more ener¬ 
getic time-scale would be of the order tl = LfU, where L 
and U are the plume diameter and velocity at the vent, which 
gives the characteristic turnover time of the largest eddies 
in a turbulent plume (e.g., Zhou et al. 2001| l. On the other 
hand, the smallest time-scale (largest Sts) can be defined by 
the Kolmogorov similarity law by ~ T^Re^ , where the 
macroscopic Reynolds number is defined, at a first instance, 
by Re^ = ULjv, v being the kinematic viscosity of the gas 
phase alone. Eor numerical models, it is also useful intro¬ 
ducing the Large-Eddy Simulation (EES) time-scale t^, with 
respect to the resolved scales related to the numerical grid 
resolution, size of the explicit filter, and discretization accu- 


racy (|Lesieur et al. 2005[ Gamier et al.| 

2009 Balachandarl 

land Eaton| 2010 Cerminara et al. 

2015 

1 . At EES scale 


Sts is not as large as at the Kolmogorov scale, thus the de¬ 
coupling between particles and the carrier fluid is mitigated 
by the EES filtering operation. We found that Sts < 0.2 for 
EES of volcanic ash finer than 1 mm. 

The model presented here is conceived for resolving di¬ 
lute suspensions, namely mixtures of gases and particles with 
volumetric concentration ^ = Cs 5,10“^. We here use the 


definition of dilute suspension by Elg hobasTu] \\99\ 1994|l 
and Balachandar| pOO^ , corresponding to regimes in which 
particle-particle collisions can be disregarded. This can also 
be justified by analyzing the time scale of particle-particle 
collisions. In the dilute regime, in which we can assume 
an equilibrium Maxwell distribution of particle velocities, 
the mean free path of solid particles is given by (|Gidaspow 
|T994l l; 

1 d. 


^'’■'’“6y2es 


(9) 


Consequently, particle-particle collisions are relatively in¬ 
frequent (Ap.p ^0.1 m ^ ds), so that we can neglect, as 
a first approximation, particle-particle collisions and con¬ 
sider the particulate fluid as pressure-less, inviscid and non- 
conductive. 
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In volcanic plumes the particle volumetric concentration 
can exceed of one order of magnitude the threshold Cs ~ 10“^ 
only near the vent (see, e.g.. Sparks et al. 1997 ?). However, 
the region of the plume where the dilute suspension require¬ 
ment is not fulfilled remains small with respect the size of 
the entire plume, weakly influencing its global dynamics. In¬ 
deed, as we will show in Sec. air entrainment and particle 
fallout induce a rapid decrease of the volumetric concentra¬ 
tion. On the contrary, the mass fraction of the solid parti¬ 
cles can not be considered small, because particles are heavy: 
es*Ps = Ps — Pg- Thus, particles inertia will be considered in 
the present model: in other words, we will consider the two 
way coupling between dispersed particles and the carrier gas 
phase. 


Summarizing, our multiphase model focuses and carefully 
takes advantage of the hypotheses characterizing the follow¬ 
ing regimes: heavy particles (ps /Pg ^ 1) in dilute suspension 
(Cs < 10“^) with dynamical length scales much larger than 
the particles diameter (point-particle approach) and relative 
Reynolds number smaller than 10^. 


2.1 Eulerian-Eulerian multiphase flow model 


When the Stokes number is smaller than one, and the num¬ 
ber of particles is very large, it is convenient to use an Eule¬ 
rian approach, where the carrier and the dispersed phase are 
modeled as interpenetrating continua, and their dynamics is 


described by the laws of fluid mechanics (Balachandar and 
|Eaton|[20T0| l. 

Here we want to model a polydisperse mixture of 
i S [1,2,...,/] = 7 gaseous phases and j G [1,2,..., J] = (J 
solid phases. From now on, we will use the subscript {■)j 
instead of (•)s for the jth solid phase. Solid phases represent 
the discretization of a virtually continuous grain-size distri¬ 
bution into discrete bins, as usually done in volcanological 
studies (cf. |Cioni et ak] |2003[ |Neri et ak] |2003| l. Another 
possible approach is the method of moments, in which the 
evolution of the moments of the grain size distribution is de¬ 
scribed. This has recently been applied in volcanology to in¬ 
tegral plume models by de’ Michieli Vitturi et al. ( 2015| l. In 
the present work we opted for the classical discretization of 
the grain size distribution (cf. |Neri et al.||200^ . In ( |Cermi-| 
nara[ |20I5b|l we analyze the Eulerian-Eulerian model under 


the barotropic regime to show the existence of weak solutions 
of the corresponding partial differential equations problem. 


' dtPi + y ■{piUi)=Q, 

^tPj “f V ■ {Pj'^j ) ~ ; j ^ 3: 

(PgUg) + V • (pgUg 0 Mg) -f Vp = 

= y-T + pg-J2fj; 

dt {pjUj) + V- {pjUj 0 Mj) = 

~ Pj9 f j 1 J G 3; 

dt {Pghg) + V • (pghgUg) -f V • (q - T • Mg) = 

= dtp-dtipgKg) - V • (pg/fgMg)-^ 

-|-pg(q• Mg) — {uj ■ f j + Qj)', 

je3 

^dtipjhj) -f V • {pjhjUj) = Qj +Sjhj , j G 3; 


with the following constitutive equations (g is the gravita¬ 
tional acceleration): 

- Given the mass fractions of the gaseous (solid) 
phases and pm the bulk density of the mixture, 
the bulk density of the gas phase is Pg = ^jPi = 
ThiViPm, while the mass fraction of the solid phases 
Ps — ^^ 3 Pj ~ Pm- Consequently, pj^ = pg -t-ps. 

- The volumetric concentration of the ith(jth) phase is 
given by ei=pilpi. 

- Perfect gas: p = J2jPiRiTg, with Ri the gas constant of 
the ith gas phase. This law can be simplified by nothing 
that Cs 1, thus Ci ~ 1 and pi ~ pi (cf. [Suzuki] |2005| ). 
Anyway, in this work we use the complete version of 
the perfect gas law. It can be written in convenient form 
for a poly-disperse mixture as: 

1 ^y^yj I 

Tto P 


- Newtonian gas stress tensor: 

T = 2p(Tg)(sym(VMg)-^V-Mgl), (12) 

where p{T) — is the gas dynamic viscosity 

and Pi is that of the ith gas component. 

- Enthalpy per unit of mass of the gas (solid) phase: 

= with Cp^i (Cj) the 

specific heat at constant pressure of the ith (jth) phase. 

- The Fourier law for the heat transfer in the gas: 
q = —k^VT, where fcg = 'Y^^eiki and ki is the conduc¬ 
tivity of the ith gas component. 


In the regime described above, the Eulerian-Eulerian equa¬ 
tions for a mixture of a gas and a solid dispersed phase are 


(jFeireislj |2004[ [Marble] [1970] [Neri et al.[ [2003 1 [Gidaspow] 
1 1994[[Gamier et al.[[2009[[Berselli et al.|[2015t ?): 


- Qj and fj refer to Qt and fd of the jth solid phase; 
Sj is the source or sink term (when needed) of the jth 
phase. Ki = |Mi|^/2 is the kinetic energy per unit of 
mass of the ith gas phase (Kj for the jth solid phase). 
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2.2 Equilibrium-Eulerian model 


while for the phases we have: 


In the limit St^ <C 1, the drag terms fj and the thermal ex¬ 
change terms Qj can be calculated by knowing Mg and Tg, 
and the Eulerian-Eulerian model can be largely simplified 
by considering the dusty-gas (also known as pseudo-gas) ap¬ 
proximation ( [Marble 19701. A refinement of this approx¬ 
imation (valid if Stj < 0.2), has been developed by Maxey 
( |TW| l, as a first-order approximation of the Lagrangian par¬ 
ticle momentum balance (see Eq. ([TOlid): 

dtUj+Uj-S/Uj =—{ug — Uj)+g. (13) 


By using the Stokes law and a perturbation method, and by 
defining a = D^Mg (with D* = ^ -f m • V ), we obtain a cor¬ 
rection to particle velocity up to first order 


Uj = Ug + Wj — Tj{a + Wj ■\7ug) + 0{Tj) (14) 


leading to the so-called equilibrium-Eulerian model de¬ 
veloped by Eerry and Balachandar ([2001 |l and Balachan- 


jdar and Eatcm (20101 for incompressible multiphase flows. 
It is worth noting that at the zeroth order we recover 
Uj =Ug + Wj, where Wj is the settling velocity defined in 
Eq. (0. 


To write the compressible version of that model, we define 
the relative jth particle velocity held Vj so that Uj —Ug + Vj. 
Recalling the definitions of the mass fraction and the mixture 
density given above, we define: 


u, = -^yjVj (15) 

u^ — Ug — u^ (16) 

T^ = '^{yjVj®Vj)-u^®u^, (17) 


By summing up the gas momentum equation with the solid 
momentum equations in Eq. ([TOll, we thus obtain: 


(PiritTm) ^ ' (Pm^m ^ “b Pm^r) — 

= -\/p + 'V ■T + p^g + '^SjUj. (18) 
i63 

This momentum balance equation is equivalent to the com¬ 
pressible Navier-Stokes equation with the substitution Mg — 
Mm and the addition of the term V • (pmTTr) which takes into 
account the first order effects of particle decoupling on mo¬ 
mentum (two-way coupling). We keep this term because of 
the presence of the settling velocity Wj in Vj which is at the 
leading order. 

Moving to the mass conservation, summing up over i and 
j the continuity equations in ( fTO) !, we obtain the continuity 
equation for the mixture: 

9tPm + V-(pmMm) = 

J63 


9t(pm2/i) + V- (pmMgy*) = 0, ieo (20) 

(PmPj ) ~b V • [pm(Mg)?/j] = , j G. 3 • (71) 


It is worth noting that the mixture density follows the clas¬ 
sical continuity equation with velocity held Mm. We refer to 
Mm as the mixture velocity held. 

As pointed out in Eq. ([^ and below, in our physical regime 
the thermal Stokes time is of the same order of magnitude of 
the kinematic one. However, this regime has been thoroughly 


analyzed in the incompressible case by Eerry and Balachan- 


|dar| p005| l, demonstrating that the error made by assum¬ 
ing thermal equilibrium is at least one order of magnitude 
smaller than that on the momentum equation (at equal Stokes 
number), thus justifying the limit Tj -^Tg = T as done for 
the thermal equation in the dusty gas model. 

By summing up all enthalpy equations in ( [TOl i, and 
by defining = 

we obtain; 


(Pm^m) H” ^ ' [Pm^m ("^m H” ’^c)] — 

— ^tP (Pm-^m) ^ ’ [Pm^m ("^m H” 

+ y ■{T-u^-q)+p^{g-u^) + '^Sj{hj+Kj). (22) 

ie3 


The terms 


Vc=U,+ 
Vk=Ui- + 






Km 


(23) 

(24) 


take into account the combined effect of the kinematic de¬ 
coupling and the difference between the specific heat (vc) 
and kinetic energy (vx) of the mixture and of the jth specie. 

Summarizing, the compressible equilibrium-Eulerian 
model is: 


^tPm “b ^ ' (Pm^tTm) — ^ ^ ? 

je3 

“b V • — 0 , t € 5, 

(PniPj ) “b V • {^PijiUj yj ) = Sj , j G 3 : 

dt (PmMm) + V • (pmMtn ® Mm -b Pm^r) = 

= -Vp + V -T + p^g + y^SjUj; 

363 


(Pm^m) H” ^ ’ [Pm^m('^m ’^c)] — 

= dtp - dt (Pmi^m) - V • [pmi^m('Wm + )] + 

■ {T■ Ug - q) + p^{g ■ u^) +'^Sj{hj + Kj). 

J63 


(25) 

The first equation is redundant, because it is contained in the 
second and third set of continuity equations. Note that we 
have not used the explicit form of Vj for deriving Eqs. ([251), 
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which therefore can be used for any multiphase flow model 
with I phases moving with velocity Ug and temperature 
T, and J phases each moving with velocity Uj =Ug + Vj 
and temperature T. However, in what follows we will use 
Eq. ( |T4l i when referring to the compressible Equilibrium- 
Eulerian model. 

It is also worth noting that in the Navier-Stokes equa¬ 
tions it is critical to accurately take into account the non¬ 
linear terms contained by the conservative derivative + 
V • (V’w) because they are the origin of the major difficul¬ 
ties in turbulence modeling. A large advantage of the dusty 
gas and equilibrium-Eulerian models is that in both models 
the the most relevant part of the drag fj) and heat ex¬ 
change (^g Qj) terms have been absorbed into the conserva¬ 
tive derivatives for the mixture. This fact allows the numeri¬ 
cal solver to implicitly and accurately solve the particles con¬ 
tribution on mixture momentum and energy (two-way cou¬ 
pling), using the same numerical techniques developed in 
Computational Eluid Dynamics for the Navier-Stokes equa¬ 
tions. The dusty gas and Equilibrium-Eulerian models are 
best suited for solving multiphase system in which the parti¬ 
cles are strongly coupled with the carrier fluid and the bulk 
density of the particles is not negligible with respect to that 
of the fluid. 

The equilibrium-Eulerian model thus reduces to a set of 
mass, momentum, and energy balance equations for the gas- 
particle mixture plus one equation for the mass transport of 
the particulate phase. In this respect, it is similar to the dusty- 
gas equations, to which it reduces for Ts = 0. With respect to 
the dusty-gas model, here we solve for the mixture velocity 
ttm, which is slightly different from the carrier gas velocity 
Mg. Moreover, here kinematic decoupling is taken into ac¬ 
count by moving the I gas phases and the J solid phases 
with different velocity fields, respectively Mg and Uj. Thus, 
we are accounting for the imperfect coupling of the particles 
to the gas flow, leading to preferential concentration and set¬ 
tling phenomena (the vector Vj includes a convective and a 
gravity accelerations terms). 

The equilibrium-Eulerian method becomes even more ef¬ 
ficient (relative to the standard Eulerian) for the polydis- 
perse case (J > 1). Eor each species of particle tracked, 
the standard Eulerian method requires four scalar fields, the 
fast method require one. Eurthermore, the computation of 
the correction to Vj needs only to be done for one particle 
species. The correction has the form —Tja, so once the term 
a is computed, velocities for all species of particles may be 
obtained simply by scaling the correction factor based on the 
species’ response times Tj. To be more precise, the standard 
Eulerian method needs /-I-5J-I-4 scalar partial differential 
equations, while the equilibrium-Eulerian model needs just 
I + J + A, i.e. 4J equations less. 


2.3 Sub-grid scale models 


The spectrum of the density, velocity and temperature fluctu¬ 
ations of turbulent flows at high Reynolds number typically 
span over many orders of magnitude. In the cases where the 
turbulent spectrum extend beyond the numerical grid resolu¬ 
tion, it is necessary to model the effects of the high-frequency 
fluctuations (those that cannot be resolved by the numerical 
grid) on the resolved flow. This leads to the so-called Large- 
Eddy Simulation (LES) techniques, in which a low-pass Alter 
is applied to the model equations to Alter out the small scales 
of the solution. In the incompressible case the theory is well- 
developed (see Berselli et al.[ 2005 [ Sagaut 2006| l, but LES 
for compressible flows is still a new and open research held. 
In our case, we apply a spatial Alter, denoted by an overbar 
(i5 is the Alter scale): 


i/;= f G{x — x';6)^{x')dx'. 

Ja 


(26) 


Some example of LES Alters G{x;S) used in compressible 
turbulence are reviewed in [Gamier et al.| ( |2009| ). In com¬ 
pressible turbulence it is also useful to introduce the so-called 
Eavre Alter: 




Puli’ 

Pm 


(27) 


Eirst, we apply this Alter to the Equilibrium-Eulerian 
model fundamental equation ([T4|) modifled as follows: 


Uj =Ug+ Wj -\- 

- (5t Mm -f Mtn • VM™ + (MJr -f MJj ) • VMm) -f O(rj) (28) 

moving all the new second order terms into 0{tJ), using 
dtUj -f Uj • Vt/j = 0 and deflning: 

w, = -'^yjWj. (29) 

3 

Multiplying the new expression for Uj by and Eavre- 
flltering, at the flrst order we obtain: 


Pm'^j — pm (^g 4“ j ) 4” 

T? (Pm'^m) 4“ V • (pm'^m ^ ’^m) 4“ Pm('^T 4“ Wj ) • VMm) 4“ 

— (30) 


where we have used Tj = Tj and consequently Wj = Wj be¬ 
cause the Stokes time changes only at the large scale and it 
can be considered constant at the Alter scale. Moreover, we 
have deflned the subgrid-scale Reynolds stress tensor: 

B = (31) 


As discussed and tested in Shotorban and Balachandaij 
(20071, the subgrid terms can be considered 0{Tj) and ne¬ 
glected when multiplied by flrst order terms. Another form 
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of Eq. ( |30l l can be recovered by noting that at the leading 
order {tm ~ Eg — - 11 ),: 


iij — iig+Wj — Tj (dtUg + iig-Vug + Wj • VEg +V -B/pn,). 

(32) 

We recall here the Boussinesq eddy viscosity hypothesis: 


a 


(33) 


where the deviatoric part of the subgrid stress tensor can be 
modeled with an eddy viscosity pt times the rate-of-shear 
tensor Sm = sym(VEni) — |V • Eml- The first term on the 
right hand side of Eq. ( |3^ is the isotropic part of the subgrid- 
scale tensor, proportional to the subgrid-scale kinetic energy 
Ki. While in incompressible turbulence the latter term is ab¬ 
sorbed into the pressure, it must be modeled for compressible 
flows (cf Moin et aT] ( |19911 l and |Yoshizaw a|( [T9861 l). |Ducros| 
et al. ( 1995|l showed another way to treat this term by absorb¬ 


ing it into a new macro-pres sure and macro-temperature (cf. 
also [Lesieur et al.| ( |2005| l and [Lodato et al.| ( |20()9[ l). We re¬ 


call here also the eddy diffusivity viscosity model (cf. Moin 


|et al. ( 19911 l): any scalar 'ip transported by Mm generates a 
subgrid-scale vector that can be modeled with the large eddy 
variables. We have: 


Pm (MmV’- , 


(34) 


where Pp is the subgrid-scale turbulent Prandtl number. 

We apply the Eavre filter defined in Eq. ( |27] l to Eqs. 

(for the application of the Eavre filter to the compressible 
Navier-Stokes equations cf. Gamier et al. (|2009 1 , Moin et al. 


( 1991| l and Erlebacher et al. 


(19901), obtaining: 


l9tPm + V- (pmMm) =^m; 

V • (pmMgP^) = 

{PmVj ) T V • [pmMj ijj ] = Sj V • 0 j , j G 3 ■, 

(Pm^tn) ^ ' (Pm'tJ’m ^ “t" PmTr) “b Vp — 

= V • T-f ^ Sj'iij -j-p^g — V • B 

dtiPm^m) “b V ■ [Pm(^m ~b Mpr)^m] — 

= - l9t (PmiTm ) - V • [pm (Em-b ) .Am] + 

+ V-(f •Eg-g)-bPm(g-Em) + 

+ J2s,{h,+K,)-V-Q, 


(35) 

(36) 

(37) 

(38) 




where 


— PmiVi'ttg Vi'^g) — 

= PmiyjUj — yi'iij) = 

2 

B = Pm ('tZm ® tZm '^m ® ’^m ) — 3 Pm Atl 2ptSn 

d 

Q — Pm((^m’trm ^rntZ-m) — ^ ; 

rp 


(39) 

(40) 

(41) 

(42) 

(43) 


are: the subgrid eddy diffusivity vector of the ith phase; of 
the jth phase; the subgrid-scale stress tensor; the diffusiv¬ 
ity vector of the temperature; respectively. Other approx¬ 
imations have been used to derive the former EES model: 
the viscous terms in momentum and energy equations, and 
the pressure-dilatation and conduction terms in the energy 
equations are all non-linear terms and we here treat them as 
done by Erlebacher et al. ( | 1990) 1 and Moin et al. ( 1991) 1. The 
subgrid terms corresponding to the former non-linear terms 
could be neglected so that, for example, p V • Ug ~ p V • Eg. 
In particular, this term has been neglected also in presence 
of shocks (cf. )Garnier et al.) ( )2002l l). We refer to )Vrem^ 
( 1995) 1 for an a priori and a posteriori analysis of all the ne¬ 
glected terms of the compressible Navier-Stokes equations. 
Moreover, in our model the mixture specific heat Cm and 
the mixture gas constant i?m vary in the domain because 
Ui and Uj vary. Thus, also the following approximations 
should be done, coherently with the other approximations 
used: = C^T ~ CmE and R^T ~ RmT- 

In order to close the system, terms pi. At and Pp must 
be chosen on the basis of EES models, either static or dy¬ 
namic (see Moin et )I99I) jBardina et al.))I980[ )Germano| 
et al. I99I)l. In the present model, we implemented several 


sub-grid scale (SGS) models to compute the SGS viscosity, 
kinetic energy and Prandtl number (Cerminara 2015b I. Cur¬ 
rently, the code offers the possibility of choosing between: 
1) the compressible Smagorinsky model, both static and dy- 
namic (see jPureby 1996) Yoshizawa 199^ Pope) )2000) 
)Chai and Mahesh) )2()12) )Gamier et al.) 2009) 1 ; 2) the sub¬ 
grid scale K-equation model, both static and dynamic (see 


Chacon-Rebollo and Lewandowski 

12013) Fureby )1996) 

Yoshizawa) ) 1993) )Chai and Mahesh 

2012)l; 3) the dynami- 

cal Smagorinsky model in the form b 
the WALE model, both static and dy 

y Moin et al. ()1991 1; 4) 

namic (seejNicoud and] 

Ducros)]1999)]Lodato et al.)]2009)]Piscaglia et al. 

20131. 


All through this paper, we present results obtained with 
the dynamic WALE model (see Eig.j^and the corresponding 
section for a study on the accuracy of this EES model). A 
detailed analysis of the influence of subgrid-scale models to 
simulation results is beyond the scopes of this paper and will 
be addressed in future works. 


3 Numerical solver 

The Eulerian model described in Section is solved nu¬ 
merically to obtain a time-dependent description of all in¬ 
dependent flow fields in a three-dimensional domain with 
prescribed initial and boundary conditions. We have cho¬ 
sen to adopt an open-source approach to the code develop¬ 
ment in order to guarantee control on the numerical solu¬ 
tion procedure and to share scientific knowledge. We hope 
that this will help building a wider computational volcanol¬ 
ogy community. As a platform for developing our solver, 
we have chosen the unstructured, finite volume (EV) method 
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based open source C++ library OpenFOAM (version 2.1.1). 
OpenFOAM, released under the Gnu Public License (GPL), 
has gained a vast popularity during the recent years. The 
readily existing solvers and tutorials provide a quick start 
to using the code also to inexperienced users. Thanks to a 
high level of abstraction in the programming of the source 
code, the existing solvers can be freely and easily modified 
in order to create new solvers (e.g., to solve a different set 
of equations) and/or to implement new numerical schemes. 
OpenFOAM is well integrated with advanced tools for pre¬ 
processing (including meshing) and post-processing (includ¬ 
ing visualization). The support of the OpenCFD Ltd, of 
the OpenFOAM foundation and of a wide developers and 
users community guarantees ease of implementation, main¬ 
tenance and extension, suited for satisfying the needs of both 
volcanology researchers and of potential users, e.g. in vol¬ 
cano observatories. Finally, all solvers can be run in paral¬ 
lel on distributed memory architectures, which makes Open¬ 
FOAM suited for envisaging large-scale, three-dimensional 
volcanological problems. 

The new computational model, called ASHEE (ASH 
Equilibrium Eulerian model) is documented in the 
VMSG (Volcano Modeling and Simulation Gate¬ 
way) at Istituto Nazionale di Geofisica e Vulcanologia 
(http;//vmsg.pi.ingv.it) and is made available through the 
VHub portal (https;//vhub.org). 


3.1 Finite Volume discretization strategy 


In the FV method ( [Ferziger and Peri3| 1996| l, the govern¬ 
ing partial differential equations are integrated over a com¬ 
putational cell, and the Gauss theorem is applied to convert 
the volume integrals into surface integrals, involving surface 
fluxes. Reconstruction of scalar and vector fields (which are 
defined in the cell centroid) on the cell interface is a key step 
in the FV method, controlling both the accuracy and the sta¬ 
bility properties of the numerical method. 

OpenFOAM implements a wide choice of discretization 
schemes. In all our test cases, the temporal discretization is 


based on the second-order Crank-Nicolson scheme (Ferziger 


and Peric I996| l, with a blending factor of 0.5 (0 meaning a 
first-order Euler scheme, 1 a second-order, bounded implicit 
scheme) and an adaptive time stepping based on the max¬ 
imum initial residual of the previous time step (|Kay et aT 


2010| l, and on a threshold that depends on the Courant num¬ 
ber (C < 0.2). All advection terms of the model are treated 
implicitly to enforce stability. Diffusion terms are also dis¬ 
cretized implicit in time, with the exception of those repre¬ 
senting subgrid turbulence. The pressure and gravity terms 
in the momentum equations and the continuity equations are 
solved explicitly. However, as discussed below, the PISO 
(Pressure Implicit with Splitting of Operators, |Issa| ( |1986| l) 
solution procedure based on a pressure correction algorithm 
makes such a coupling implicit. Similarly, the pressure ad¬ 
vection terms in the enthalpy equation and the relative veloc¬ 


ity Vj are made implicit when the PIMPLE (mixed SIMPLE 


and PISO algorithm, Ferziger and Peric ( |I996| l) procedure 
is adopted. The same PIMPLE scheme is applied treating 
all source terms and the additional terms deriving from the 
equilibrium-Eulerian expansion. 

In all described test cases, the spatial gradients are dis¬ 
cretized by adopting an unlimited centered linear scheme 
(which is second-order accurate and has low numerical dif¬ 
fusion - Ferziger and Peric 1996|l. Analogously, implicit 


advective fluxes at the control volume interfaces are recon¬ 
structed by using a centered linear interpolation scheme (also 
second order accurate). The only exception is for pres¬ 
sure fluxes in the pressure correction equation, for which 
we adopt a TVD (Total Variation Diminishing) limited lin¬ 
ear scheme (in the subsonic regimes) to enforce stability and 


non-oscillatory behavior of the solution. We refer to Jasak 
( 1996[ ) for a complete description of the discretization strat 
egy adopted in OpenFOAM. 


3.2 Solution procedure 


Instead of solving the set of algebraic equations deriving 
from the discretization procedure as a whole, most of the 
existing solvers in OpenFOAM are based on a segregated 
solution strategy, in which partial differential equations are 
solved sequentially and their coupling is resolved by iterating 
the solution procedure. In particular, for Eulerian fluid equa¬ 
tions, momentum and continuity equation (coupled through 
the pressure gradient term and the gas equation of state) are 
solved by adopting the PISO algorithm. The PISO algorithm 
consists of one predictor step, where an intermediate veloc¬ 
ity field is solved using pressure from the previous time-step, 
and of a number of PISO corrector steps, where interme¬ 
diate and final velocity and pressure fields are obtained it¬ 
eratively. The number of corrector steps used affects the 
solution accuracy and usually at least two steps are used. 
Additionally, coupling of the energy (or enthalpy) equation 
can be achieved in OpenFOAM through additional PIMPLE 
iterations (which derives from the SIMPLE algorithm by 
|Patankar||I980| l. For each transport equation, the linearized 
system deriving from the implicit treatment of the advection- 
diffusion terms is solved by using the PbiCG solver (Precon¬ 
ditioned bi-Conjugate Gradient solver for asymmetric matri¬ 
ces) and the PCG (Preconditioned Conjugate Gradient solver 
for symmetric matrices), respectively, preconditioned by a 
Diagonal Incomplete Lower Upper decomposition (DILU) 
and a Diagonal Incomplete Cholesky (DIC) decomposition. 
The segregated system is iteratively solved until a global tol¬ 
erance threshold epiMPLE is achieved. In our simulations, we 
typically use cpimple < 10~^ for this threshold. 

The numerical solution algorithm is designed as follows: 

1. Solve the (explicit) continuity equation ( [35] l for mixture 
density (predictor stage: uses fluxes from previous 
iteration). 
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2. Solve the (implicit) transport equation for all gaseous 

and particulate mass fractions: i = 1 ,...,/ and 

Vj, j = 

3. Solve the (semi-implicit) momentum equation to obtain 
Mm (predictor stage: uses the pressure field from previ¬ 
ous iteration). 

4. Solve the (semi-implicit) enthalpy equation to update 
the temperature field T and the compressibility Pm/p 
(pressure from previous iteration). 

5. Solve the (implicit) pressure equation to update the 
pressure (uses predicted values of fluxes). 

6. Correct density, velocity and fluxes with the new pres¬ 
sure held (keeping T and Pm/p fixed). 

7. Iterate from [^evaluating the continuity error as the dif¬ 
ference between the kinematic and thermodynamic cal¬ 
culation of the density (PISO loop). 

8. Compute explicit terms (transport coefficients and de¬ 
coupling). 

9. Evaluate the numerical error cpimple and iterate fromj^ 
if prescribed (PIMPLE loop). 

10. Compute EES subgrid terms. 

With respect to the standard solvers implemented in 
OpenEOAM (v2.1.1) for compressible fluid flows (e.g. 
sonicFoam or rhoPimpleFoam), the main modification 
required are the following: 

1. The mixture density and velocity replaces the fluid ones. 

2. A new scalar transport equation is introduced for the 
mass fraction of each particulate and gas species. 


4 Verification and validation study 


A wide set of numerical tests has been performed to assess 
the adequacy of the ASHEE model for the intended vol- 
canological application and the reliability of the numerical 
solution method. Validation tests are focused on the dynam¬ 


ics of gas (Section 4.1 1 and multiphase (Section [4.2| l turbu¬ 
lence and on the mixing properties of buoyant plumes (Sec¬ 
tion |4^. Compressibility likely exerts a controlling role to 


the near-vent dynamics during explosive eruptions (e.g., Car- 
cano et akl 2013| l. Although this is not the focus of this 
work, we briefly discuss in Section 4.4 the performance of 
the model on a standard one-dimensional shock wave numer¬ 
ical test. 


4.1 Compressible decaying homogeneous and isotropic 
turbulence 


Turbulence is a key process controlling the dynamics of vol¬ 
canic plumes since it controls the rate of mixing and air en¬ 
trainment. To assess the capability of the developed model to 
resolve turbulence (which requires low numerical diffusion 
and controlled numerical error [Geurts and Erohli^ |2002[ ), 
we have tested the numerical algorithm against different con¬ 
figurations of decaying homogeneous and isotropic turbu¬ 
lence (DHIT). 


In this configuration, the flow is initialized in a domain 
n which is a box with side L = 2tt with periodic boundary 


conditions. As described in Lesieur et al. (200511; Honein and| 

Moin (200411; 

Liao et al.|(|20091l; Pirozzoli and Grasso (2004]l; 

Blaisdell et a 

. (19911, we chose the initial velocity field so 


that its energy spectrum is 


£(fc) 



(44) 


3. The equations of state are modified as described in 
Eq.([TT). 

4. Eirst-order terms from the equilibrium-Eulerian model 
are added in the mass, momentum and enthalpy equa¬ 
tions. 

5. Equations are added to compute flow acceleration and 
velocity disequilibrium. 

6. Gravity terms and ambient fluid stratification are added. 

7. New SGS models are implemented. 


Concerning point! 


Eerry et al. (2003 


intj^ it i; 
)03]i7the 


is worth remarking that, accordingly to 
first-order term in Tj in Eq.([T^ must 
be limited to avoid the divergence of preferential concentra¬ 
tion in a turbulent flow field (and to keep the effective Stokes 
number below 0.2). In other word, we impose at each time 
step that \Tj{a + Wj • VMg)| < 0.2|Mg-f |. We tested the 


effect of this limiter on preferential concentration in Sec. 4.2 
below. 


with peak initially in k = ko and so that the initial kinetic 
energy and enstrophy are: 


1 

Ko = J £(A)dfc=(45) 

poo r 

:Ko= / k^E{k)dk=-u^^^kQ. (46) 

Jo o 


As reviewed by Pope ( 2000[ ), the Taylor microscale can be 
written as a function of the dissipation e = 2i/3-C: 


2 _ _ 5K 

/\x ^ - 


e 7f ’ 

thus in our configuration, the initial Taylor micro scale is: 

l5Ko 2 


(47) 


At,o — 


7f( 


(48) 


0 


^0 


We have chosen the non-dimensionalization keeping the root 
mean square of the magnitude of velocity fluctuations (m') 
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equal to Ui-ms: 

We also chose to make the system dimensionless by fixing 
Pm,o = To = Pr = so that the ideal gas law becomes: 

P = /9mi?mT’ = T?m, (50) 


pOO 

a/ u' ■ u'dx = 2 / £(fc)dfc. 

Jo 


(49) 


and the initial Mach number of the mixture based on the ve¬ 
locity fluctuations reads: 


Ma^ms 



I2Kqp„ 

ImP 


= Mrms(7mP)’ 


(51) 


This means that Marms can be modified keeping fixed Unns 
and modifying p. Following |Honein and Moin] ( |2004| l, we 
define the eddy turnover time: 


Te = 


^rms 


(52) 


The initial compressibility ratio Co is defined as the ratio 
between the kinetic energy and its compressible component 
Kc: 


T^c,o 


Ko 


1 

2(2^)3iTo 



■u'^dx. 


(53) 


Here, u'^ is the compressible part of the velocity fluctuations, 
so that V -u' — V -u'^ and V A= 0. 

The last parameter, i.e. the dynamical viscosity, can be 
given both by fixing the Reynolds number based on At,o or 
ko: 


Rex = 
Refco = 


PmMrmsAT,0 

\/3p 


Pm^rms 

kop. 


(54) 

(55) 


Ma^s = 0.2, Co = 0, = 2Ko = 0.056, ko =4, At = 

0.5, Te ~ 3.6596, p = 5.885846 * 10-^ Re^ ~ 116, Refc„ ~ 
100. Thus a grid with N = 2563 cells gives kj^ax — 127 
and kmaxVK = 27r, big enough to have a DNS. The simula¬ 
tion has been performed on 1024 cores on the Fermi Blue 
Gene/Q infrastructure at Italian CINECA super-computing 
center (http://www.cineca.it), on which about 5 h are needed 
to complete the highest-resolution runs (2563 of cells) up to 
time //re = 5.465 (about 3500 time-steps). The average re¬ 
quired total CPU time on 1024 Fermi cores is about 1-3 mil¬ 
lions of cells per second, with the variability associated with 
the number of solid phases described by the model. This 
value is confirmed in all benchmark cases presented in this 
paper. 

Fig. [3 reports the parallel efficiency on both the Fermi 
and the PLX (a Linux cluster based on Intel Xeon esa- and 
quad-core processors at 2.4 GHz) machines at CINECA. The 
ASHEE code efficiency is very good (above 0.9) up to 512 
cores (i.e., up to about 30000 cells per core), but it is overall 
satisfactory for 1024 cores, with efficiency larger than 0.8 on 
PLX and slightly lower (about 0.7) on Fermi probably due to 
the limited level of cache optimization and input/output seal- 
ability (|Culpo 2011| l. The code was run also on 2048 cores 
on Fermi with parallel efficiency of 0.45 (Dagna 20131. 



PLX 


It is useful to define the maximum resolved wavenumber on 
the selected A^-cells grid and the Kolmogorov length scale 
based on Re^g. They are, respectively: 


(?-; 

. 27r N 

(56) 

' L N-1' 

TT^ -3 

-Re,;. 


(57) 


Fig. 1: ASHEE parallel efficiency on Fermi and PLX super¬ 
computers at CINECA (www.cineca.it). 

Fig. 1^ shows an isosurface of the second invariant of the 
velocity gradient, defined as: 

Q„ = i(Tr(VM)2-(VM.VM)), (58) 


In order to have a DNS, the smallest spatial scale 5 should 
be chosen in order to have k^axilK > 2 (jPirozzoli and Grassoj 
[2004^ . 

We compare the DNS of compressible decaying homoge¬ 
neous and isotropic turbulence with a reference, well tested 
numerical solver for Direct Numerical Simulations of com¬ 
pressible turbulence by Pirozzoli and Gras^pO()4| l; |Bernar-| 
|dini and Pirozzoli (20091. For this comparison we fix the 
following initial parameters: p = = 1, 7ni = 1.4, Pr = 1, 


The so called Q-criterion ( [Gamier et al.||200^ allows indeed 
the identification of coherent vortices inside a three dimen¬ 
sional velocity field. 

In Fig. we present a comparison of the energy spec¬ 
trum E{k) obtained with the ASHEE model and the model 
by [Bernardini and Pirozzoli] ( |2009| l after approximatively 1 
eddy turnover time; the norm of the difference between 
the two spectra is 4.0 * 10“^. This validates the accuracy of 
our numerical code in the single-phase and shock-free case. 
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Fig. 2: Isosurface at — 19 Hz^ and f/rg ~ 2.2, represent¬ 
ing zones with coherent vortices. 



Fig. 3: Comparison of a DNS executed with the eight or¬ 
der scheme by Prrozzoli and Grasso ( |2004| l and our code im¬ 
plemented using the C++ libraries of OpenFOAM at t/r^ = 
1.093. The norm between the two spectra is 4.0* 10“^. 
The main parameters are Re^ — 116, Manns = 0.2. 


in Fig. 18 and 19 by [Gamier et al. (1999 1 . 

Fig. 4b shows the kinetic energy spectrum at t/T^ = 


0,1.093,5.465. We notice that the energy spectrum widens 
from the initial condition until its tail reach k ~ fcmax — 127. 
Then system becomes to dissipate and the maximum of the 
energy spectrum decreases. The largest scales tend to lose 
energy slower than the other scales and the spectrum widens 
also in the larger scale direction. 

Fig. I^presents the evolution of K (total turbulent kinetic 
energy), TC (enstrophy). At (Taylor microscale). We notice 
that the total kinetic energy decreases monotonically and at 
t ~ 5.5re just ~ 15% of its initial value is conserved. On the 
other hand, enstrophy increases until it reaches a maximum 
at 1.5 < f/Te < 2. It then starts to decrease monotonically. 
This behavior is related to the two different stages we have 
highlighted in the analysis of the energy spectrum evolution. 
In the first stage, viscous effects are negligible and enstrophy 
increases due to vortex stretching. During the second stage, 
viscous diffusion starts to have an important role and dis¬ 
torted dissipative structures are created (Gamier et al. 1999| ). 
Also the Taylor microscale reflects this behavior, reaching a 
minimum at the end of the first stage and increasing mono¬ 
tonically during the second stage of the evolution. It is a 
characteristic of the magnitude of the velocity gradients in 
the inertial range: by comparing it with 6 we can have an 
idea of the broadness of the range of wave numbers where 
the flow is dissipative. In this DNS, we have At — 10.2^ at 
t ~ 5.5Te. 

In Fig. 1^ we show the evolution of the Kolmogorov time 
scale tk during the evolution of the decaying turbulence. 

We finally compare in Fig.j^the DNS described with sim¬ 
ulations at lower resolution with N = 32^ and N = 64^ cells. 
In this case, it is expected that the spectra diverge from the 
DNS, unless an appropriate subgrid model is introduced to 
simulate the effects of the unresolved to the resolved scales. 


Several subgrid models have been tested (Cerminara 2015b I, 


both static and dynamic. Fig. [^presents the resulting spec¬ 
trum using the dynamic WAFF model ( jNicoud and Ducro^ 
1999| Fodato et al.j 20091. In this figure, we notice how the 
dynamic WAFF model works pretty well for both the 32^ 
and 64^ FES, avoiding the smallest scales to accumulate un¬ 
physical energy. 


4.2 Two-phase isotropic turbulence 


Fig. 0 shows the evolution of several integral param¬ 
eters describing the dynamics of the decaying homoge¬ 


neous and isotropic turbulence. Fig. 4a displays the den¬ 
sity fluctuations = sj{{p-{p)ny)u, the density con¬ 
trast Pmax/Pmin and the standard measure of compressibil¬ 
ity C = (|V • Mp)o/(| which takes value between 0 

(incompressible flow) and 1 (potential flow) ( jBoffetta et al.j 
[2007^ . All the quantities showed in Fig. [^depend on the ini¬ 
tial Mach number end compressibility. For the case shown, 
Marms = 0.2 and we obtain very similar result to those ported 


In this section we test the capability of our numerical code to 
correctly describe the decoupling between solid and gaseous 
phases when St^ < 0.2 and to explore its behavior when the 
equilibrium-Eulerian hypothesis Stj < 0.2 is not fulfilled so 
that a limiter to the relative velocity Ug — Uj is applied. 

To this aim, we performed a numerical simulation of ho¬ 
mogeneous and isotropic turbulence with a gas phase initial¬ 
ized with the same initial and geometric conditions described 


in Sec. 4.1 We added to that configuration 5 solid particle 
classes (j = 2-+6) chosen in such a way that Stj G [0.03,1], 
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Fig. 4: Evolution of some dynamical quantities in DHIT with Re^ — 116 and Marms = 0.2 at f/rg = 5.465. 



Fig. 5: Energy spectrum £(fc) at f/xg = 5.465 obtained with 
different spatial resolutions and with/without subgrid scale 
LES model. 


homogeneously distributed and with zero relative velocity; 
Vj(x,0) = 0. From Fig. 4d we see that, during turbulence 


decay, approximately tk €[0.6,1.2]. Therefore, for a given 


Tj 

Stmax ^Tj/O.e 

^3 (pj — 10^) 

0.60 

1.0 

2.521*10“^ 

0.30 

0.5 

1.783*10“^ 

0.15 

0.25 

1.261*10“^ 

0.075 

0.125 

8.914* 10-^ 

0.0375 

0.0625 

6.303* 10-^ 


Table 1; Stokes time, maximum Stokes number and diameter 
of the solid particles inserted in the turbulent box. 


particle class with tj fixed, during the time interval tlr^G 
[0,5.5] we have St^ax/Stmin — 2. In Tab.[^we report the main 
properties of the particles inserted in the turbulent box. To 
evaluate the Stokes time here we used tj = pjd'^/{18^) be¬ 
cause in absence of settling Re^ < 1 when St^ < 1 (jBalachan- 
|dar| |2009| l. We set the material density of all the particles 
to pj = 10^. In order to have a small contribution of the 
particle phases to the fluid dynamics - one way coupling - 
here we set the solid particles mass fraction to a small value, 
Pj = 0.002, so that j/g = 0.99. 

In Fig. 1^ we show a slice of the turbulent box at f/rg ~ 
2.2. Panel a) displays the solid mass fraction, highlighting 
the preferential concentration and clustering of particles in 
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(a) Mass fraction (b) Acceleration 


Fig. 6 : Slice of the turbulent box at t/r^ ~ 2.2. The two 
panels represent respectively a logaritmic color map of 2/3 
(Stmax = 0.5) andof |ag|. 


response to the development of the acceleration field (panel 
b) associated with turbulent eddies. 

As described in Maxey (|1987[); Rani and Balachandar 


( |2003| l, a good measure for the degree of preferential con¬ 
centration in incompressible flows is the weighted average 
on the particle mass fraction of the quantity (|I)p — |§p), 
where § is the vorticity tensor, i.e. the skew symmetric part 
of the gas velocity gradient and T) is its symmetrical part. 
For compressible flows, we choose to consider 


(T), = ((|B|2-|§|2-|Tr(2))p)),- = 

(y-(y)o))o 

{yj)n 


(59) 


This is a good measure because (use integration by parts. 
Gauss theorem and Eq. ([T^ with Wj = 0): 


\ l^m 

= -r,((|D| 2 -|§| 2 _|Tr( 2 ))| 2 ))^. (60) 

Moreover, it is worth noting that (CP)j vanishes in absence of 
preferential concentration. By dimensional analysis, prefer¬ 
ential concentration is expected to behave as: 



(T)j cx 




DNS 

FES, 


( 61 ) 


because it must be proportional to Tj and have a dimension 
of [s“^]. As described by [Pope 2001], the typical time-scale 
corresponding to an eddy length-scale ^ in the inertial sub¬ 
range, can be evaluated by means of the Kolmogorov’s the¬ 
ory as: 


T^=TX 



2 

3 


(62) 



Fig. 7: Evolution of the degree of preferential concen¬ 
tration with Stj (FES) or St^ (DNS). We obtain a good 
agreement between equilibrium-Eulerian FES/DNS and Fa- 


grangian DNS simulations. The fit for the data by Rani and 


Balachandar 


(20031 is found in Eq. 


where the Taylor microscale At is defined by Eq. 47 Since 
the time based on the Taylor microscale is defined as 


TX = 


•\/3At 


(63) 


we can evaluate the typical time at the smallest resolved EES 
scale ^ knowing the kinetic energy K{t) and AT(f): 


Tj(f) = . 


2K{t) 




(64) 


In Fig.j^we show the time-evolution of the degree of pref¬ 
erential concentration as a function of the Stokes number for 
both DNS with 256^ cells and the EES with 32^ cells. There, 
we multiply (CP)j by in order to make it dimensionless 
and to plot on the same graph all the different particles at 
different times together. 

At f = 0 the preferential concentration is zero for all Stokes 
number. Then, preferential concentration of each particle 
class increases up to a maximum value and then it decreases 
because of the decaying of the turbulent energy. The max¬ 
imum degree of preferential concentration is reached by 
each particle class when tk is minimum (at f/xe ~ 1.7, cf. 


Fig. 4d I. Then, (CP) j decreases and merges with the curve 


relative to the next particle class at the final simulation time, 
when Tk is about twice its minimum. Note that the expected 
behavior of Eq. ( |M] l is reproduced for St^ <0.2 and in par¬ 
ticular we And: 




1.52Stjr],. ■ 
1.52Stjr-- 


DNS 

EES. 


(65) 


Moreover, by comparing our results with the Eulerian- 


Lagrangian simulation described in Rani and Balachandar 
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( 2003| l, we note that our limiter for the preferential concen¬ 
tration when St > 0.2 is well behaving. 

For the sake of completeness, we found that the best ht in 


the range St < 2.5 for the data found by Rani and Balachan- 
|dar| ( |200^ is: 


(J>)j~1.52* 


St, 


-2 


1 + 3.1 * Stj + 3.8 * St^ 


2 'K 


( 66 ) 


with root mean square of residuals 8.5 * 10“^. 

For what concerns the 32^ LES simulation. Fig. shows 
that the Stokes number of each particle class in the LES case 
is much smaller than its DNS counterpart. Accordingly with 
[Balachandar and Eaton |(|2OT0ll, we have 


St? = St,(^ 


(67) 


conhrming that the equilibrium-Eulerian model widens its 
applicability under the LES approximation. We also notice 
that the presented LES is able to reproduce the expected de¬ 
gree of preferential concentration with a satisfactory level of 
accuracy when St < 0.2. In particular, the LES slightly over¬ 
estimates preferential concentration and the time needed to 
reach the equilibrium and to “forget” the particle initial con¬ 
dition. 

4.3 Turbulent forced plume 

As a second benchmark, we discuss high-resolution, three- 
dimensional numerical simulation of a forced gas plume, 
produced by the injection of a gas flow from a circular inlet 
into a stable atmospheric environment at lower temperature 
(and higher density). Such an experiment allows to test the 
numerical model behavior against some of the fundamental 
processes controlling volcanic plumes, namely density varia¬ 
tions, non-isotropic turbulence, mixing, air entrainment, and 
thermal exchange. Our study is mainly aimed at assessing 
the capability of the numerical model to describe the time- 
average behavior of a turbulent plume and to reproduce the 
magnitude of large-scale fluctuations and large-eddy struc¬ 
tures. We will mainly refer to laboratory experiments by 
[George et al. (1977 1 an djShabbir and George ( [1994 1 and nu¬ 
merical simulations by [Zhou et ~ pOOlj l for a quantitative 
assessment of model results. 

Numerical simulations describe a vertical round forced 
plume with heated air as the injection fluid. The plume axis 
is aligned with the gravity vector and is subjected to a pos¬ 
itive buoyancy force. The heat source diameter 25o is 6.35 
[cm], the exit vertical velocity on the axis uq is 0.98 [m/s], 
the inflow temperature Tq is 568 [K] and the ambient air tem¬ 
perature Ta is 300 [K]. The corresponding Reynolds number 
is 1273, based on the inflow mean velocity, viscosity and di¬ 
ameter. Air properties at inlet are Cp = 1004.5 [J/(K kg)]; 
77=287 [J/(Kkg)];^ = 3x 10-5 [Pa s]. 


As discussed by Zhou et al. ( 2001| l the development of the 
turbulent plume regime is quite sensitive to the inlet condi¬ 
tions: we therefore tested the model by adding a periodic per¬ 
turbation and a non-homogeneous inlet profile to anticipate 
the symmetry breaking, and the transition from a laminar to 
a turbulent flow. The radial profile of vertical velocity has the 
form: 


Uo{r) = 


1 -tanhI — 

4Sr \bo r 


( 68 ) 


where 6r is the thickness of the turbulent boundary layer at 
the plume inlet, that we have set at Sr = O.lfco- A periodical 
forcing and a random perturbation of intensity 0.05[/o (f) has 
been superimposed to mimic a turbulent inlet. 

The resulting average mass, momentum and buoyancy flux 
are qq = 2.03 x IQ-^ [kg s”^], toq = 1-62 x 10-5[kg m s”^], 
/o = 1.81x10-3 [kgs-i]. 

The computational grid is composed of 360 x 180 x 180 
uniformly spaced cells (deformed near the bottom plane to 
conform to the circular inlet) in a box of size 12.8 x 6.4 x 6.4 
diameters. In particular, the inlet is discretized with 400 
cells. The adaptive time step was set to keep the Courant 


number less than 0.2. Based on estimates by Plourde et al. 
( |2008| ), the selected mesh refinement is coarser than the re¬ 
quired grid to fully resolve turbulent scales in a DNS (which 
would require about 720 x 360 x 360 cells). Nonetheless, this 
mesh is resolved enough to avoid the use of a subgrid-scale 
model. This can be verified by analyzing the energy spec¬ 
tra of fluctuations on the plume axis and at the plume outer 
edges. In Fig. [^we show the energy spectra of temperature 
and pressure as a function of the non-dimensional frequency: 
the Strouhal number Str = / * 26o /uq (/ is the frequency in 
[Hz]). We recover a result similar to [Plourde et al.j ( [2008| l, 
where the inertial-convective regime with the decay —5/3 
and the inertial-diffusive regime with the steeper decay —3 
are observable ( [Lis t| [T982l l. 

Model results describe the establishment of the turbulent 
plume through the development of fluid-dynamic instabili¬ 
ties near the vent (puffing is clearly recognized as a toroidal 
vortex in Fig. |^). The breaking of large-eddies progres¬ 
sively leads to the onset of the developed turbulence regime, 
which is responsible of the mixing with the surrounding am¬ 
bient air, radial spreading of the plume and decrease of the 
plume average temperature and velocity. Figure displays 
the spatial distribution of gas temperature. Mixing becomes 
to be effective above a distance of about 4 diameters. Fig¬ 
ure 1^ displays the distribution of the vorticity, represented 
by values of the invariant (Eq. [58] ). The figure clearly 
identifies the toroidal vortex associated to the first instabil¬ 
ity mode (puffing, dominant at such Reynolds numbers). We 
have observed the other instability modes (helical and mean¬ 
dering, [Lesieur et al.[ [2005 [ l only by increasing the forcing 
intensity (not shown). 


Experimental observations by George et al. (1977 1 and 


Shabbir and George (19941 reveal that the behavior of forced 
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Fig. 8; Temperature (solid) and pressure (dashed) fluctuations energy spectra; a) at a point along the plume axis 
(0, 0, 0.5715) [m]; b) at a point along the plume outer edge (0, 0.06858, 0.5715) [m]. The slopes Str”®'^^ and Str“^ are repre¬ 
sented with a thick solid and dashed line respectively. 
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Fig. 9; Three-dimensional numerical simulation of a forced gas plume at f = 10s. a) Isosurface of temperature T = 305 [K], 
colored with the magnitude of velocity, and the temperature distribution on two orthogonal slices passing across the inlet center, 
b) Isosurface of Qu = 100 [s“^] colored with the value of the velocity magnitude, and its distribution across two vertical slices 
passing through the inlet center 
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Y-Axis 


Y-Axis 


Fig. 10; Two-dimensional slice and streamlines of the velocity field; a) time-averaged velocity field; b) instantaneous velocity 
field at f = 10 s. The mean velocity field outside the plume is approximatively horizontal while in the plume it is approximately 
vertical. The region where the mean velocity field change direction is the region where the entrainment of air by the plume 
occurs. 


plumes far enough from the inlet can be well described by 
integral one-dimensional plume models ( |Morton et al.||195^ 
|Morton[ |1959| l provided that an adequate empirical entrain¬ 
ment coefficient is used. In the buoyant plume regime at this 
Reynolds number George et al. ( 1977| l obtained an entrain¬ 
ment coefficient of 0.153. 

To compare numerical result with experimental observa¬ 
tions and one-dimensional average plume models, we have 
time-averaged the numerical results between 4 and 10 s 
(when the turbulent regime was fully developed) and com¬ 
puted the vertical mass Q(z), momentum M{z) and buoy¬ 
ancy F{z) fluxes as a function of the height. To perform this 
operation, we define the time averaging operation (“) and the 
radial domain 

I (^/tracer (^) > 0.01 ?/tracer,o) ^ (^) ^ 0)} , (69) 


hypothesis is tested in Fig. 10 1 , where it can be verified that 
the time-averaged streamlines inside the plume are parallel 
to the axis (Fig.[T^ shows the instantaneous streamlines and 
velocity magnitude field). 


The plume fluxes are evaluated as follows (cf. George 


et al.l |1977t |Shabbir and Georgej |1994| |Kaminski et al. 

2 ^ ; 


- mass flux Q(-z) = / puzdxdy 

Jfi 

- momentum flux M{z) = / puldxdy 

Jq 

- buoyancy flux F{z)= / Uz {pa — p) dxdy 

Jo. 

where pa = Pa{z) is the atmospheric density. From these 
quantities it is possible to retrieve the main plume parameters 


where (x,y,z) = a; are the spatial coordinates, i/tracer is the 
mass fraction held of a tracer injected from the vent with ini¬ 
tial mass fraction ytracer,o Uz is the axial component of the 
velocity held. We use this definition for n(z) for coherence 
with integral plume models, where the mean velocity held 
is assumed to have the same direction of the plume axis (cf. 
[Morton et al.||1956t|Woods[ri988[|Cerminara||2015a|b| l. This 


- plume radius b{z) = 

- plume density /3(z) = p^ 

- plume temperature Tfj(z) =Ta 

- plume velocity U{z) = ^ 
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- entrainment coefficient k{z) = —— 


277pet Ub 


where (•)' is the derivative along the plume axis and is the 
atmospheric temperature profile. 

Figure [TT] displays the average plume radius and velocity. 
As previously reported by |Fannel0p and Webber] ( [2003| l and 
[Plourde et al.|P00^ , the plume radius initially shrinks due to 
the sudden increase of velocity due to buoyancy (at z = 0.1 
m). Above, turbulent mixing becomes to be effective and 
increases the plume radius while decreasing the average ve¬ 
locity. The upper inset in Fig. [^represents the values of the 
vertical mass q = QjQo, momentum m = M/Mq and buoy¬ 
ancy / = F/Fq, normalized with the inlet values. All vari¬ 
ables have the expected trends and, in particular, the buoy¬ 
ancy flux is constant (as expected for weak ambient stratifica¬ 
tion) whereas q and m monotonically increase and attain the 
theoretical asymptotic trends shown also in Fig. Indeed, 
|Fannel0p and Webber ( |2003| l have shown that an integral 
plume model for non-Boussinesq regimes (i.e., large density 
contrasts) in the approximation of weak ambient stratifica¬ 
tion and adopting the Ricou and Spalding ( |1961| l formula¬ 
tion for the entrainment coefficient, has a first integral such 
that q^ is proportional to at all elevations. Figure 12 


demonstrates that this relationship is well reproduced by our 


numerical simulations, as also observed in DNS by Plourde 
|erar| ( |2008] l. 


The lower inset in Fig. 11 shows the computed entrain¬ 
ment coefficient, which is very close to the value found in 
experiments ( [George et al.||1977 Shabbir and Georg^|1994| l 
and numerical simulations ( Zhou et al.[ 2001) 1 of an analo- 
gous forced plume. We found a value around 0.14 in the 
buoyant plume region (6.4 < z/26o < 16). 

The analysis of radial profiles led to a similar conclusions: 
in Fig.[T^ we show the evolution of the radial profiles for the 
mean vertical velocity field. In this figure, we also report the 
plume radius as evaluated from Gaussian fits of these profiles 
on horizontal slices: 


= C/fitexp - 


"fit 


(70) 


The slope of the function has been evaluated in the 

region 6.4 < z/2b(j < 16, to obtain = 0.142± 0.001 to 
be compared with the result of George et al. ( 1977|: z = 

0.135±0.010. 


Finally, figure 14 reports the time-average values of the 
vertical velocity and temperature along the plume axis. As 
observed in laboratory experiments, velocity is slightly in¬ 
creasing and temperature is almost constant up to above 4 
inlet diameters, before the full development of the turbu¬ 
lence. When the turbulent regime is established, the decay 
of the velocity and temperature follows the trends predicted 
by the one-dimensional theory and observed in experiments. 
The insets displays the average value of the vertical velocity 
and temperature fluctuations along the axis. Coherently with 


experimental results ( [George et al. 19771, velocity fluctua¬ 
tions reach their maximum value and a stationary trend (cor¬ 
responding to about the 30% of the mean value) at a lower 
height (about 3 inlet diameters) with respect to temperature 
fluctuations (which reach a stationary value about the 40% 
above 4 inlet diameters). 

4.4 Transonic and supersonic flows 

Although not essential in the present application, the ability 
of solving transonic and supersonic regimes is also required 
for the full-scale simulation of volcanic processes. We here 
test the behavior of the ASHEE code in presence of shocks in 
the classical Sod’s shock tube test case ( [Sod[| 19781 ) describ¬ 
ing the expansion of a compressible, single-phase gas hav¬ 
ing adiabatic index 7 = 1.4. At f = 0 the domain of length 
10 m is subdivided in two symmetric subsets. In the first 
subset (spatial coordinate x < 0) we set u = 0, p = 10® Pa, 
T = 348.432 K, so that p=l. In the second subset (x > 0), 
we set M = 0, p = 10^ Pa, T = 278.746 K, so that p = 0.125 
kg/m®. We indicate with c = 374.348 m/s the speed of sound 
of the gas in the x < 0 part of the domain. We impose zero 
gradient boundary conditions (9a;(-) = 0 ) for all the variables 


u, p, T. As described in Sod (1978), a reference analytic 
solution exists for this problem. 

In Fig. [T^ we show the density profile obtained with 
the ASHEE model after 0.007 s of simulation. We per¬ 
formed two simulations at different resolution. The first has 
100 cells and it is compared with the OpenFOAM solver 
rhoCentralFoam with a second order semi-discrete, non 
staggered central scheme of [Kurganov et ah ( 2001[) for the 
fluxes, and a total variation diminishing limiter ([\^m^j^ 


]^97j) for the interpolation. We refer to [Greenshields et al. 
(2010) for a presentation of rhoCentralFoam and of the 


Sod’s shock tube test case. The inset of Fig. 15 is the simu¬ 
lation with an higher resolution (1000 cells). In this figure, 
we notice that the code performs satisfactorily both at low 
and high resolution. It is capable to capture the shocks pretty 
well, with a diffusion that is comparable with that obtained 
with rhoCentralFoam, a solver conceived for simulating 
shocks. 


5 3D simulation of a turbulent volcanic plume 

Numerical simulations of volcanic plumes were conducted in 
the framework of the lAVCEI (International Association of 
Volcanology and Geochemistry of the Earth Interior) plume 
model intercomparison initiative (Costa et al. 2015| ), con¬ 
sisting in performing a set of simulations using a standard 
set of input parameters so that independent results could be 
meaningfully compared and evaluated, discuss different ap¬ 
proaches, and identify crucial issues of state of the art mod¬ 
els. We here discuss three-dimensional numerical simulation 
of a weak volcanic plume in a stratified, calm atmosphere. 
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plume radius [m] 

0 0.1 0.2 0.3 0.4 0.5 



Fig. 11: Time-averaged plume radius and velocity. The insets display the non-dimensional mass, momentum and buoyancy 
fluxes (top) and the time-averaged entrainment coefficient. The line in the entrainment panel is a constant fit, from which 
results A: = 0.141 ±0.001. 



Fig. 12: Finear regression between and for the 

plume simulation with azimuthal forcing. 



r/2b„ 

Fig. 13: Radial profiles of the time-averaged velocity field at 
various height. The scale for these profiles is indicated by 
the up-down arrow on the left in the panel. The thick solid 
line is the plume radius evaluated from the mass, momentum 
and buoyancy fluxes, while the thick dashed line is the plume 
radius evaluated from Gaussian fits of horizontal profiles. 
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z/2bQ 


Fig. 14: Centerline time-average axial velocity, and temperature profiles with azimuthal forcing. Inset) centerline correlations 
of fluctuating velocity and temperature. 


whose input data were set assuming parameters and meteo¬ 
rological conditions similar to those of the 26 January 2011 
Shinmoe-dake eruption ( [Suzuki and Koyaguchi|[20T3l ). Ini¬ 
tial conditions and injection parameters are reported in Table 

m 


Parameter 

Value 

Vent elevation 

1500 m 

Vent diameter 

54 m 

Mass eruption rate 

1.5 X 10® kg/s 

Exit velocity 

135 m/s 

Exit temperature 

1273 K 

Exit water fraction 

3 wt.% 

Mixture density at vent 

4.85 kg/m® 


Table 2: Vent conditions for the weak volcanic plume simu¬ 
lation. 

The particle size distribution is composed of two individ¬ 
ual classes of pyroclasts in equal weight proportion repre¬ 
senting, respectively, hne (diameter d = 0.0625 mm; den¬ 


sity p = 2700 kg/m^, volume fraction e = 0.00086821) and 
coarse ash (diameter d = 1.0000 mm; density p = 2200 
kg/m^, volume fraction e = 0.00106553). With respect to the 
laboratory benchmark case of Section 4.3 volcanic plumes 
are characterized by non-Boussinesq regimes at the vent and 
buoyancy reversal (with the initial mixture density about 4 
times larger than the atmospheric one) and by a stratihed at¬ 
mosphere (Fig. [Tbl l. However, the most relevant difference is 
due to the significant temperature contrast (900 K) and to the 
presence of a high particle content which may strongly affect 
the mixing properties of the plume. 


The Stokes number of the solid particles is, in general, 
a complex function of time and space, since the turbulent 
flow is characterized by a wide spectrum of relevant time 
and length scales. Generally, the Stokes number is associ¬ 
ated with the most energetic turbulent eddy scale which, for 
laboratory plumes, has a typical turnover time of the order 
of tl Str|^ « 0.12 s, where Dy and are the plume 
diameter and velocity at the vent, respectively, and Str is 
the Strouhal number, of the order Str = 0.3 (jZhou et al.j 
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Fig. 15; The Sod’s shock tube density after 0.007 s (here c = 
374.348 m/s). Here we compare the analytic solution (solid) 
with two simulations performed with ASHEE model (dashed 
line) and the OpenFOAM rhoCentralFoam solver (cir¬ 
cles). The resolution is 100 cells, while in the inset it is 
reported the solution obtained with the ASHEE model with 
resolution of 1000 cells. 



density (kg/m3j pressure IhPa] temperature [K] 


Fig. 16: Atmospheric profiles as provided by the Japan Me¬ 
teorological Agency’s Non-Hydrostatic Model (| Hashimot^ 
|et al.|[20T^ for Shinmoe-dake volcano at 00 JST of 27 Jan¬ 
uary 2011. 


2001^ . Based on this time scale, and computing the particle 
relaxation time from Eq. the Stokes number for the two 
adopted particle classes is about Skoarse ~ 5 and St^ne ~ 0.2, 
so we expect to see non-equilibrium phenomena for both 
particles classes, with more evident effects on the coarsest 
phase. However, the Stokes number, as an average value in 
all the plume is not as high as calculated above. Indeed, by 
using Eq. ( |64l i as reference time for the turbulent dynamics, 
we obtain Skoarse ~ 0.6 and Stfjne ~ 0.03. It is worth recall¬ 
ing here that the equilibrium-Eulerian approach is accurate 
and advantageous for particles having St < 0.2 and that, in 
our model, we numerically limit the acceleration field in or¬ 
der to keep the turbulent non-equilibrium within this limit, 
as explained in Sect. |^and tested in Sect. |4.2| Fig. [7] The 
averaged value of this limit - measuring the importance of 
the decoupling limiter for this simulation - is approximately 
0 . 6 . 


The computational domain is cylindrical and is extended 
4836o X 7656o in the radial and vertical directions {bo be¬ 
ing the vent radius). The numerical grid is non-uniform 
and non-orthogonal. The discretization of the vent is rep¬ 
resented in Fig. ( fTTal i. For the highest resolution run, the 
cell size increases from a minimum grid size Ar = 2bo /32 
with no radial grading factor in the region where the plume 
is expected to develop (Fig. 17b| l, whose initial radius is 
equal to 2.5&o and increases linearly with an angle 9 such 
that tan6* = 0.147, slightly larger than tan6* = 0.12 predicted 
by the Morton’s plume theory with entrainment /c = 0.1 
( [Ishimine 2006| l. Outside this region, a radial grading fac¬ 
tor of 1.0446 is applied. Along z, 2048 cells are utilized. 
The minimum vertical cell size is Az = 26o/32, and a grad¬ 
ing factor of 1.00187 is imposed. The azimuthal resolution 
is constant and equal to (5.625 degrees). The resulting 
total number of cells is 10,747,904. This numerical mesh 
guarantees accuracy of the results: the solution procedure 
utilizes 2 PISO and 2 PIMPLE loops to achieve an absolute 
residual cpimple = 10“^ (see Sec.|^. 


Simulation of 720 s of eruption required about 490,000 
time steps (imposing a CEL constrain of 0.2, resulting in an 
average time-step dt « 1.5 ms, with a maximum velocity at 
the vent of about 150 m/s) for a total ran-time of about 25 
days on 1024 cores on the Fermi architecture at CINECA 
(meaning about 2.25 millions of cells per second, consis¬ 
tently with estimates of Sec.|^. 


Figure shows the development of the volcanic plume 
at / = 400 s. Because of the atmospheric stratification, the 
plume reaches a neutral buoyancy condition at about 10 km 
above the vent (i.e., 11.5 km above the sea level, still within 
the troposphere). Due to its inertia, the plume reaches its 
maximum plume height Umax ~ 12 km, higher than the neu¬ 
tral buoyancy level, before spreading radially to form the so- 
called volcanic umbrella. The two orthogonal sections high¬ 
light the different spatial distribution of the volumetric frac¬ 
tion of fine (right) and coarse (left) ash particles, due to the 
different coupling regime with the gas phase. Coarse par¬ 
ticles has indeed a larger settling velocity Wg = which 
causes a more intense proximal fallout from the plume mar¬ 
gins and a reduced transport by the umbrella. This is high¬ 
lighted by the plot of the streamlines of the mixture velocity 
along a vertical section (Fig. 19 1 , showing that the plume 
updraft is surrounded by a shell of settling coarse particles, 
which also inhibit air entrainment while promoting particle 
re-entrainment into the plume. 


Besides settling, the large inertia of the coarse ash is re¬ 
sponsible for the kinematic decoupling, leading to preferen¬ 
tial concentration and clustering of particles at the margins 
of turbulent eddies. To illustrate this phenomenon, in a non- 
homogeneous flow, the instantaneous preferential concentra¬ 
tion is computed as the (normalized) ratio between the jth 
particle concentration and the concentration of a tracer (in 
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(a) (b) 

Fig. 17: Zoom of the computational grid used for volcanic plume simulations. 



Fig. 18: Three-dimensional numerical simulation of a weak volcanic plume, 400 s after the beginning of the injection (inlet 
conditions as in Table|^. Isosurface and vertical sections of the fine (light white) and coarse (light sand) ash volume fractions. 
The two-dimensional plots represent the distribution of the volume concentration of coarse (left) and fine (right) particles 
across vertical orthogonal slices crossing the plume axis. 


our case, water vapor), i.e.. 


where the 0 subscript corresponds to the value at the vent. 


e.= 


Uj 2/tracer,0 
Vj.O ytracer 


(71) 


Fig. 20 shows the distribution of Gj for the coarsest par¬ 
ticles at f = 400 s. The color scale is logarithmic and sym¬ 
metric with respect to 1, which corresponds to the nil pref- 
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Fig. 19: Vertical section of the instantaneous value of the mixture velocity modulus (in logarithmic scale) at f = 400 s and 
velocity streamlines. 


erential concentration. For Cj < 1, the mixture is relatively 
depleted of particles (green to blue scale); for e, > 1, parti¬ 
cles are clustered (green to red scale), with mass fraction up 
to 5 times larger and 20 times smaller than the value it would 
have in absence of preferential concentration. This behavior 
is expected to affect the mixing and entrainment process. It 
is also worth remarking that the more uniform red area be¬ 
yond the plume margins corresponds to the region of settling 
particles below the umbrella region. On the other hand, the 
top of the plume is relatively depleted of coarse particles. 
The corresponding Figure]^ for fine particles confirms that 
these are tightly coupled to the gas phase and almost behave 
as tracers (value of Cgne is everywhere around 1). These con¬ 
clusions are coherent with the a-priori estimate of St^ we 
gave at the beginning of this section, based on the Taylor mi¬ 
croscale time (Eq. (|64|i). 


Finally, we present the results obtained by averaging the 
volcanic plume flow field over time (in a time-window 
[300-720] s where the plume has reached statistically sta¬ 
tionary conditions) and over the azimuthal angle, in order 
to allow comparison with one-dimensional integral models 
(e.g.. Woods 1988| l and discuss the effect of numerical reso¬ 
lution. The averaging procedure is the generalization of that 
explained in Sect. I^to the multiphase case (see Cerminara 


|2015a| ). The form of the results presented are similar to those 
obtained in Fig. [^for the laboratory plume test case. 

Figure [22| presents the results of the averaging procedure 


for three multiphase flow simulations at different resolution 
(panels a-c). In particular, panel a) has the highest resolu¬ 
tion (minimum radial cell size Ar = 2&o/32 with 2bo equal 
to the inlet diameter); panel b) has Ar = 26o/16; panel c) 
has Ar = 2bo/8. In panel d) we present results at the lowest- 
resolution obtained by imposing the full kinematic equilib¬ 
rium between gas and particles, i.e., by adopting the dusty- 
gas model (Marble 1970| l. 


Results demonstrate that the numerical model is quite ro¬ 
bust and accurate so that even low-resolution simulations are 
able to capture the main features of the volcanic plume de¬ 
velopment. However, the maximum plume height systemat¬ 
ically decreases from 12100 m (a), to 11300 m (b) to 11000 
m (c) when we decrease the resolution. Analogously, the 
Neutral Buoyancy Level (NBL) decreases from 7800 m (a) 
to 7200 m (b) to 7100 m (c). Although the lowest resolu¬ 
tion run seems to underestimate the maximum plume height 
and the plume radius by about 10%, the average velocity pro¬ 
file (and the vertical profiles of q, m and /) is consistent in 
the three runs, showing a jet-plume transition at about 2000 
m above the vent, also corresponding to the transition to a 
super-buoyancy region ( Woods [ 2010| l. The computed en¬ 
trainment coefficient is also consistent and relatively inde¬ 
pendent on the grid resolution and shows a different behav¬ 
ior with respect to the laboratory case, associated with the 
effect of the density contrast. In this case, a maximum value 
of about fc ~ 0.1 is obtained in the buoyant plume region be- 
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Fig. 20: Distribution of Ccoarse (Eq. 71 1 for the coarsest particles across a vertical section at f = 400 s. 
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Fig. 21: Distribution of Cfjne (Eq. 71 1 for the finest particles across a vertical section at f = 400 s. 
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tween 2 and 5 km above the vent. 

Interestingly, we find that in three-dimensional simula¬ 
tions the entrainment decreases near the NBL and it become 
negative above that level. This happens because the mass exit 
from the plume region defined in Eq. ( [69] l moving from it to 
the umbrella cloud. In this way, the mass flow q of the plume 
decreases above the NBL and a stationary solution can be 
achieved. This is not the case in integral plume models with 
positive entrainment coefficient, where the maximum plume 
height is reached as a singularity point with divergent mass 
flow and plume radius (cf. |Mortonl|1959||Woods|['l988| l. We 
plan to study this behavior more thoroughly in future studies. 

The dusty-gas model shows a significantly different be¬ 
havior, with a larger plume radius, a slightly higher entrain¬ 
ment coefficient and a more marked jet-plume transition with 
no further acceleration (without a super buoyancy transition). 
The plume height is slightly lower than the non-equilibrium 
case at the same resolution having maximum plume height 
and neutral buoyancy level of 9900 m and 6100 m, respec¬ 
tively. Numerical simulations thus suggest that the effects 
of non-equilibrium gas-particle processes (preferential con¬ 
centration and settling) on air entrainment and mixing are 
non-negligible. These effects are certainly overlooked in the 
volcanological literature and will be studied more thoroughly 
in future studies, by applying the present model to other re¬ 
alistic volcanological case studies. 

6 Conclusions 

We have developed a new, equilibrium-Eulerian model to nu¬ 
merically simulate compressible turbulent gas-particle flows. 
The model is suited to simulate relatively dilute mixtures 
(particle volume concentration e < 10“^) and particles with 
Stokes number St < 0.2. It is appropriate to describe the dy¬ 
namics of volcanic ash plumes, with kinematic decoupling 
between the gas and the particles, assumed in thermal equi¬ 
librium. 

We have tested the model against controlled experiments 
to assess the reliability of the physical and numerical for¬ 
mulation and the adequacy of the model to simulate the 
main controlling phenomena in volcanic turbulent plumes, 
and in particular: 1) multiphase turbulence (including prefer¬ 
ential concentration and density effects); 2) buoyancy and 
compressibility effects; 3) stratification and density non¬ 
homogeneity. 

The model reproduces the main features of volcanic 
plumes, namely: 1) buoyancy reversal and jet-plume tran¬ 
sition; 2) plume maximum height and spreading of the um¬ 
brella above the neutral buoyancy level; 3) turbulent mixing 
and air entrainment; 4) clustering of particles; 5) proximal 
fallout and re-entrainment of particles in the plume. Results 
demonstrate that the compressible equilibrium-Eulerian ap¬ 
proach adopted in the ASHEE model is suited to simulate 
the three-dimensional dynamics of volcanic plumes, being 


able to correctly reproduce the non-equilibrium behavior of 
gas-particle mixtures with a limited computational cost. 

Einally, the adopted open-source computational infras¬ 
tructure, based on OpenEOAM, will make the model easily 
portable and usable and will ease the maintenance and im¬ 
plementation of new modules, making ASHEE suitable for 
collaborative research in different volcanological contexts. 
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Fig. 22: Time-averaged plume radius and velocity. The insets display the non-dimensional mass, momentum and buoyancy 
fluxes (top) and the time-averaged entrainment coefficient. Panels a-c) ASHEE model at different resolutions; panel d) dusty- 
gas model. 
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